#Table of contents {.tabset .tabset-fade .tabset-pills}

##Sample distribution As each locus follows binomial distribution, the genetic drift can be modelled \(\frac{\sqrt{pq}}{2n_e}\), in which \(n_e\) is the effective population size.

f=0.5
pop=c(50, 200, 500, 1000)
gn=100
cohort=100
G=matrix(0, cohort, gn)
layout(matrix(1:4, 2, 2))
for(p in 1:length(pop)) {
  plot(main=paste("Population size", pop[p]), x=NULL, y=NULL, xlim=c(1, gn), ylim=c(0, 1), xlab="Generation", ylab="Frequency", bty='n')
  for(i in 1:cohort) {
    G[i,1] = 0.5
    for(j in 2:gn) {
      G[i,j]=mean(rbinom(pop[p], 2, G[i,j-1]))/2
    }
    lines(G[i,], col=sample(1:10, 1))
  }
}

Subgroups

Three pop

M=c(10000)
fst1=0.05
fst2=0.01
for(m in 1:length(M)) {
  P=runif(M[m], 0.2, 0.8)
  Sz=c(100, 100, 100, 100, 100, 100)
  Fq=matrix(0, length(Sz), M[m])
  Fq[1,]=rbeta(M[m], P*(1-fst1)/fst1, (1-P)*(1-fst1)/fst1)

  P3=rbeta(M[m], P*(1-fst1)/fst1, (1-P)*(1-fst1)/fst1)

  P2=(0.5*Fq[1,]+0.5*P3)

  Fq[3,]=rbeta(M[m], P2*(1-fst2)/fst2, (1-P2)*(1-fst2)/fst2)
  Fq[4,]=rbeta(M[m], P2*(1-fst2)/fst2, (1-P2)*(1-fst2)/fst2)
  
  Fq[2,]=rbeta(M[m], P3*(1-fst2)/fst2, (1-P3)*(1-fst2)/fst2)
  Fq[5,]=rbeta(M[m], P3*(1-fst2)/fst2, (1-P3)*(1-fst2)/fst2)
  Fq[6,]=rbeta(M[m], P3*(1-fst2)/fst2, (1-P3)*(1-fst2)/fst2)

  G=matrix(0, sum(Sz), M[m])
  for(i in 1:Sz[1]) {
    G[i,] = rbinom(M[m], 2, Fq[1,])
  }
  
  for(i in (Sz[1]+1):(Sz[1]+Sz[2])) {
    G[i,] = rbinom(M[m], 2, Fq[2,])
  }
  
  for(i in (Sz[1]+Sz[2]+1):sum(Sz)) {
    G[i,] = rbinom(M[m], 2, Fq[3,])
  }

  for(i in (Sz[1]+Sz[2]+Sz[3]+1):sum(Sz)) {
    G[i,] = rbinom(M[m], 2, Fq[4,])
  }

  for(i in (Sz[1]+Sz[2]+Sz[3]+Sz[4]+1):sum(Sz)) {
    G[i,] = rbinom(M[m], 2, Fq[5,])
  }

  for(i in (Sz[1]+Sz[2]+Sz[3]+Sz[4]+Sz[5]+1):sum(Sz)) {
    G[i,] = rbinom(M[m], 2, Fq[6,])
  }

  Gs=apply(G, 2, scale)
  GG=Gs %*% t(Gs)
  EigenG=eigen(GG)
  layout(matrix(1:2, 1, 2))
  barplot(EigenG$values[1:20])
  plot(EigenG$vectors[,1], EigenG$vectors[,2])
  layout(matrix(1:6, 2, 3, byrow=T))
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[1:100,1], EigenG$vectors[1:100,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[101:200,1], EigenG$vectors[101:200,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[201:300,1], EigenG$vectors[201:300,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[301:400,1], EigenG$vectors[301:400,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[401:500,1], EigenG$vectors[401:500,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
  plot(xlim=range(EigenG$vectors[,1]), ylim=range(EigenG$vectors[,2]), main=paste(M[m], "markers"), EigenG$vectors[501:600,1], EigenG$vectors[501:600,2], bty='n', xlab="PC 1", ylab="PC 2", pch=16)
}

Wishart distribution simulation

Generate Wishart distribution

## [1]   10   10 1000
##             [,1]       [,2]        [,3]       [,4]        [,5]        [,6]
##  [1,]  1.0000000  0.9892037  0.85009974  0.7349126  0.40302917 -0.10687984
##  [2,]  0.9892037  1.0000000  0.90912646  0.7834566  0.41542351 -0.12399023
##  [3,]  0.8500997  0.9091265  1.00000000  0.9135839  0.57037342  0.01718272
##  [4,]  0.7349126  0.7834566  0.91358394  1.0000000  0.82871062  0.38087610
##  [5,]  0.4030292  0.4154235  0.57037342  0.8287106  1.00000000  0.79783054
##  [6,] -0.1068798 -0.1239902  0.01718272  0.3808761  0.79783054  1.00000000
##  [7,] -0.5194526 -0.5528924 -0.45133366 -0.1292602  0.32855574  0.74949700
##  [8,] -0.8033035 -0.8441922 -0.76921814 -0.5102579 -0.05009135  0.45706156
##  [9,] -0.9787197 -0.9846968 -0.85143400 -0.7074475 -0.31657491  0.17689287
## [10,] -0.9771913 -0.9736050 -0.85344880 -0.7742355 -0.44638410  0.02528311
##             [,7]        [,8]       [,9]       [,10]
##  [1,] -0.5194526 -0.80330348 -0.9787197 -0.97719132
##  [2,] -0.5528924 -0.84419222 -0.9846968 -0.97360500
##  [3,] -0.4513337 -0.76921814 -0.8514340 -0.85344880
##  [4,] -0.1292602 -0.51025791 -0.7074475 -0.77423553
##  [5,]  0.3285557 -0.05009135 -0.3165749 -0.44638410
##  [6,]  0.7494970  0.45706156  0.1768929  0.02528311
##  [7,]  1.0000000  0.87737055  0.5566532  0.43394504
##  [8,]  0.8773705  1.00000000  0.8480151  0.75653301
##  [9,]  0.5566532  0.84801514  1.0000000  0.97486002
## [10,]  0.4339450  0.75653301  0.9748600  1.00000000

##Distribution of the Wishart diagonal elements \(G=\frac{1}{m}\tilde{X}^T\tilde{X}\), and the diagonal of \(G\) follow \(\frac{\chi^2_n}{n}\)

##Homo cohort

##Heterogeneous cohort

##Tracy-Widom distribution

##TW